% Matlab code to visualize the equilibrium for the calibrated model

% To run this file seperately, specify the root directly below:
% cd 'C:/Users/Jilun/Dropbox/WorldEconomy_RiskIntolerance/paper RFS final submission/replication/numerical and calibration exercises'


% This code is for the case without LSAPs. It makes two figures
% FIGURE 6: Plot the equilibrium price for a range of productivites (included in the main text)
% FIGURE B1: Graphically illustrate the equilibrium for a particular productivity shock (included in the appendix) 


clear all;
% Baseline parameters (will change later to illustrate multiplicity)
global g rho sigma PStar taub tauh k l tau tauBarBench zStar demandImpactBench elasticityBench
loadVariables_calibration; % Sets the above parameters and the baseline shocks

zHigh = 1;
% Find the threshold that eliminates the good equilibrium
zLow = zStar-0.02;
zPost = 0.97;

ph = exp(rho+g - tauh/sigma^2);

% Recall that zh=l/ph
    %zh = l/ph; 
    %zLow= zh;
    
zArr = [zLow:0.0001:zHigh];
pArr = zArr*0;
% Solve the equilibrium for a range of z's
for i=1:length(zArr)
    z0= zArr(i);
    eq{i} = solveEquilibrium(z0,0);
    pArr(i) = eq{i}.p;
    uniqueArr(i) = eq{i}.unique; % Whether the equilibrium is unique
    % When there are multiple equilibria, the code automatically picks the
    % highest equilibrium
end

% With these parameters, there is always a bad equilibrium
% Find the first z that also allows for a good equilibrium
zInd= find(uniqueArr==0,1,'first'); 

%%%%%%%%%%%% FIGURE 6
%%%%%%%%%%%% Calibrated equilibrium and its approximation without LSAPs

figure('Position',[30 30 830 430]);
plot(zArr(zInd:end),pArr(zInd:end),'LineWidth',3);
xlabel('Productivity, $z_0$','Interpreter','Latex','FontSize',20);
title('Price per productivity, $p_0$ (in the best equilibrium)','Interpreter','Latex','FontSize',20);
zLow = zArr(zInd); zHigh =zStar+0.005;
pLow = pArr(zInd); pHigh = pArr(end);
xlim([zLow zHigh]); ylim([pLow pHigh]);
set(gca,'Xtick',[0:0.005:1]);

% Plot a line and mark zStar
line([zStar zStar],[pLow pHigh]);
specText = sprintf('$z^*$');
text(zStar,pLow+0.002,specText,'Interpreter','Latex','FontSize',20);

% Plot the pre-shock approximation line
% First calculate the slope
priceImpactPreShock = demandImpactBench/elasticityBench;
zArrPreShock = [zLow:0.001:zStar];
% Approximate in log units around zStar
pArrPreShock = exp((log(zArrPreShock)-log(zStar))*priceImpactPreShock);
hold on;
plot(zArrPreShock,pArrPreShock,'b--','LineWidth',1.5);
specText = sprintf('Approximation with $z_0p_0=1$');
text(zLow,pArrPreShock(1)+0.015,specText,'Interpreter','Latex','FontSize',20);


% Now plot the post-shock approximation line
% First calculate the slope
% Plot this for z0=zPost==0.97 (roughly 1.5% decline relative to zStar)
zIndPost= find(zArr>=zPost,1,'first'); 
pPost = pArr(zIndPost); %The equilibrium value of z0P0
% Mark this equilibrium
plot(zPost,pPost,'ko-','LineWidth',3);


demandImpactPostShock = k*(l/(zPost*pPost))/(tauh/(taub-tauh) + k*(1-(l/(zPost*pPost))));
elasticityPostShock = 1/rpBench - demandImpactPostShock;
priceImpactPostShock = demandImpactPostShock/elasticityPostShock;
zArrPostShock = zArrPreShock;
% Approximate in log units around zStar
pArrPostShock = exp((log(zArrPostShock)-log(zStar))*priceImpactPostShock);
hold on;
plot(zArrPostShock,pArrPostShock,'b:','LineWidth',1.5);
specText = sprintf('Approximation with $z_0p_0=%1.2f$',zPost*pPost);
text(zPost+0.0004,pPost-0.003,specText,'Interpreter','Latex','FontSize',20);

saveas(gcf,[pwd '/figures/fig6.eps'],'epsc');


%%%%%%%%%%%% FIGURE B1
%%%%%%%%%%%% illustrate the benchmark case (z_0=1) and a productivity shock that induces a demand recession
z0_1 = 1;
z0_2 = zPost;
eq{1} = solveEquilibrium(z0_1,0);
eq{2} = solveEquilibrium(z0_2,0);

pLow = 0.9;
% Plot for pArr>pLow

figure('Position',[30 30 830 430]);
% Plot the remaining curves
lineType.actual{1} = 'b--';
lineType.required{1} = 'r--';
lineType.actual{2} = 'b';
lineType.required{2} = 'r';

for i=1:2
% Plot each equilibria respectively
hold on;
if eq{i}.rf > 0
% The program returns the actual Sharpe ratio with rf=0
% When there is a unique corner solution
% make an adjustment to the Sharpe ratio to accommodate the fact that rf>0
sharpeActualArr = eq{i}.sharpeActualArr - eq{i}.rf/sigma;
legendStr{i*2-1} = sprintf('Actual Sharpe ratio, $r_0^f>0$');
elseif eq{i}.rf == 0
sharpeActualArr = eq{i}.sharpeActualArr;
legendStr{i*2-1} = sprintf('Actual Sharpe ratio, $r_0^f=0$');
end
sharpeRequiredArr = eq{i}.sharpeRequiredArr;
pArr=eq{i}.pArr;
ind = eq{i}.ind; % Index of the equilibrium
z0 = eq{i}.z0;

h(i*2-1) = plot(eq{i}.pArr,sharpeActualArr,lineType.actual{i},'LineWidth',3);
h(i*2) = plot(pArr,sharpeRequiredArr,lineType.required{i},'LineWidth',3);
legendStr{i*2} = sprintf('Required Sharpe ratio, $z_0=%1.2f$',z0);

% Mark ph on the figure
specText = sprintf('$p^h$');
text(ph,eq{1}.sharpeRequiredArr(end)+0.01,specText,'Interpreter','Latex','FontSize',20);
% Mark the equilibrium
if i==1 % corner equilibrium with efficient price
plot(pArr(ind),sharpeRequiredArr(ind),'ko-','LineWidth',3);    
else
plot(pArr(ind),sharpeActualArr(ind),'ko-','LineWidth',3);
end

end

titleStr = sprintf('Supply shocks with constrained interest rate in a calibrated equilibrium');
title(titleStr,'Interpreter','Latex','FontSize',20);
xlabel('Normalized price, $p_0$','Interpreter','Latex','FontSize',20);

hL = legend(h([3 4 1 2]),legendStr{3},legendStr{4},legendStr{1},legendStr{2});
set(hL,'orientation','vertical','fontsize',15,'interpreter','latex');

xlim([pLow 1]);

saveas(gcf,[pwd '/figures/figB1.eps'],'epsc');

